Microbiomes and Planctomycete diversity in large-scale aquaria habitats

In commercial large-scale aquaria, controlling levels of nitrogenous compounds is essential for macrofauna health. Naturally occurring bacteria are capable of transforming toxic nitrogen species into their more benign counterparts and play important roles in maintaining aquaria health. Nitrification, the microbially-mediated transformation of ammonium and nitrite to nitrate, is a common and encouraged process for management of both commercial and home aquaria. A potentially competing microbial process that transforms ammonium and nitrite to dinitrogen gas (anaerobic ammonium oxidation [anammox]) is mediated by some bacteria within the phylum Planctomycetes. Anammox has been harnessed for nitrogen removal during wastewater treatment, as the nitrogenous end product is released into the atmosphere rather than in aqueous discharge. Whether anammox bacteria could be similarly utilized in commercial aquaria is an open question. As a first step in assessing the viability of this practice, we (i) characterized microbial communities from water and sand filtration systems for four habitats at the Tennessee Aquarium and (ii) examined the abundance and anammox potential of Planctomycetes using culture-independent approaches. 16S rRNA gene amplicon sequencing revealed distinct, yet stable, microbial communities and the presence of Planctomycetes (~1–15% of library reads) in all sampled habitats. Preliminary metagenomic analyses identified the genetic potential for multiple complete nitrogen metabolism pathways. However, no known genes diagnostic for the anammox reaction were found in this survey. To better understand the diversity of this group of bacteria in these systems, a targeted Planctomycete-specific 16S rRNA gene-based PCR approach was used. This effort recovered amplicons that share <95% 16S rRNA gene sequence identity to previously characterized Planctomycetes, suggesting novel strains within this phylum reside within aquaria.


Introduction
The importance of microbiomes to the health of their hosts and environments is undisputed [1][2][3]. Microbiomes play critical roles in ecosystem viability by contributing to nutrient and sand filtration systems were sampled in triplicate from four different tanks (designated T7, T20, T30, and T34) differing in temperature, salinity, turnover rates and resident macrofauna (Table 1). Water turnover is mediated by large mechanical filtration systems that effectively aerate each of these aquaria. Dissolved oxygen is monitored in all tanks and does not fall below 90% saturation, showing minimal variation (<3%) across depth and time. The water in these tanks is continuously cycled through the filtration systems. These systems do not undergo full water exchanges. For water samples, approximately 1 L of surface collected water was filtered through Sterivex™ cartridges with 0.22 μm membrane filters (Millipore-Sigma, USA). For each sand filtration system, approximately 30 g of sand was collected into Whirl-Paks1 bags (Nasco, USA) from drained units using sterilized spoons. Samples for metagenomic analyses were collected from the same set of tanks and also a biological denitrification system attached to tank T30 during a single sampling time point in January 2018. Immediately following collection, all samples were frozen in liquid nitrogen until transferred to the University of Tennessee where they were stored at -80˚C until DNA extractions were conducted.

16S rRNA gene-based microbial community analysis
DNA was extracted from water and sand samples using a phenol-chloroform protocol. For water samples, the initial steps of the extraction protocol were performed within Sterivex™ cartridges by plugging the outflow port with Cha-seal clay (DWK Life Sciences, USA). Reagents were added to the cartridges directly using a needle and syringe; Luer-lock caps were used to seal the inflow port. For each sand sample, 0.5 g of material was placed in a 15 ml plastic tube (Falcon). To all samples, 1.7 ml CTAB extraction buffer, 65 μL proteinase K (10mg/ml), 65 μL lysosome (10 mg/ml), and 162 μL of filter-sterilized SDS (10% in deionized water) were added. Tubes and cartridges were incubated in a rotary agitator at 65˚C for 2 h. Following incubation, 800 μL aliquots were pipetted into 1.7 ml microcentrifuge tubes and cooled to 4˚C. The remaining solution was stored at -20˚C for later additional extraction of low yield or low-quality samples. An equal volume of phenol:chloroform:isoamyl alcohol (P:C:I, 25:24:1, pH 8.0) was added to each sample and vortexed. The aqueous and organic layers were separated via 10,000 rpm centrifugation for 5 minutes at 4˚C. The aqueous layer was transferred to a new sterile tube. Two additional P:C:I extractions were repeated with the aqueous layer. To the final aqueous layer, 450 μL of 100% isopropanol was added and incubated overnight at room temperature. Tubes were centrifuged at 10,000 rpm for 20 minutes to pellet DNA. After decanting, DNA was washed with 70% ethanol, centrifuged at 10,000 rpm for 5 minutes, and dried for 15 minutes in a laminar flow hood. DNA was then suspended in 50 μL of sterile nuclease-free water at 50˚C. DNA was quantified using a Nano-drop1 ND-1000 Spectrometer (Thermo Fisher Scientific, USA) and samples were stored at -80˚C. 16S rRNA genes were amplified using the EarthMicrobiome project PCR primers 515F-Adapt and 806R-Adapt following published protocols (S1 Table [ 27]). Amplified PCR

PLOS ONE
products were~270 bp in size, allowing for overlap in paired ends reads [27]. PCR products were sequenced using the MiSeq Illumina1 platform at the University of Tennessee's Genomic Core Facility. The Mothur software package (version 1.39.5) was used to process sequences and remove low quality reads following established criteria [28]. Mothur was also used to cluster sequences into operational taxonomic units (OTUs; at >97% identity to cultured organisms) following the Schloss MiSeq SOP for 16S rRNA genes analysis [28]. The Primer-e software package (version 7) was used to interrogate the relationships between OTUs across samples. Alpha diversity (Shannon-Weiner and Simpson indices) was also calculated using Primer-e. Pairwise Wilcoxon tests were performed to identify significant differences between diversity measures of habitat samples using the R package ggpubr (version 0.2.5; https://CRAN.R-project.org/package=ggpubr). Detection of "biomarker" OTUs, diagnostic of specific aquaria, was performed using a Linear Discriminant Analysis using the LEfSe tool as part of the BioBakery meta'omics analysis environment [29]. This analysis focused on the two parameters (substrate and salinity) contributing most to clustering of communities on the Primer ordination plot. The input data was provided as read counts and an all-vs-all analysis was performed after the substrate (water vs sand) and salinity (fresh vs salt) were collapsed.

Metagenomic sequencing and analysis
For metagenome analyses, DNA was extracted from 9 samples (water and sand filters from each exhibit and water from the T30 denitrification tank) using the MO-BIO Power Soil1 DNA isolation kit following the manufacturer's guidelines (QIAGEN, Germany). DNA was prepared for metagenomic analyses using the Nextera1 XT genome kit (Illumina1, United States) following manufacturer's protocol with the minor modification of increasing the PCR cycles from 12 to 15 because of low DNA yield in some of the samples. The samples were loaded at 8 pM with 2% PhiX spike-in on an Illumina MiSeq at the University of Tennessee Genomics Core on a v3, 600 cycle flow cell, reading 275 bases paired-end. Paired-end sequencing data was imported into CLC Genomics Workbench (version 20.0.01; QIAGEN, Germany) where reads were assembled into contigs using the default settings of the De Novo Assembly tool. Coding domains were identified using the MetaGeneMark program and annotated with CLC. Individual reads were then mapped back to assembled contigs for quantification of relative abundances of individual genes. The GhostKoala program [30] within the Kyoto Encyclopedia of Genes and Genomes (KEGG; [31]) was used to automatically annotate CDS (coding sequence) regions within each metagenome and highlight specific genes within various functional pathways.

Planctomycete-specific 16S rRNA gene amplification and analysis
Due to the relatively high proportion of Planctomycete OTUs found in sand filtration samples from T30 and T34, DNA was extracted from June 2017 archived sample material using the MO-BIO Power Soil1 DNA isolation kit (QIAGEN, Germany). A Planctomycete-specific, nested 16S rRNA gene PCR approach was employed to examine the diversity of Planctomycete species in the aquaria and probe for Planctomycetes capable of specific metabolic functions (anammox) (S1 Fig). The initial PCR amplification using the Planctomycete-specific primer Pla46F [32] and the universal bacterial primer 1390R (S1 Table [ 33]) yielding products of~1.3 kb. PCR products of the appropriate size were excised from agarose gels using a QIAquick Gel Extraction kit (QIAGEN, Germany) and used for a subsequent round of PCR amplification with the degenerate primers AMXU368F and AMXU820R (S1 Table). This primer pair is diagnostic of Planctomycetes capable of performing anammox and generates an expected product size of~470 bp [34]. PCR products were cleaned using a QIAquick PCR Purification Kit (QIAGEN, Germany), directly ligated into the pCR4.0™-TOPO™ vector (Invitrogen, USA), and transformed into chemically competent E. coli cells (One Shot™ TOP10; Invitrogen, USA). Recombinant clones were selected on LB agar supplemented with ampicillin (50 μg/ml). Colony PCR was performed on isolated colonies by centrifuging 100 μl of turbid culture for 3 minutes at 13,000 rpm. The pellets were suspended in 100 μl Milli-Q water and then incubated for 10 minutes at 95˚C. Following centrifugation for 2 minutes at 10,000 rpm, the supernatant was transferred to a new, sterile, 0.5 ml tube and stored at -20˚C. Colony PCR was performed using M13F and M13R primers that recognize the cloning vector. Twenty-five colonies were screened for each sample type (T30 and T34 sand filters) and insert sizes of both~470 bp and 1.3 kb were obtained. Plasmid extractions were performed with a QIAprep1 Spin Miniprep kit (QIAGEN, Germany) for all clones and sent to the University of Tennessee's Genomics Core for Sanger sequencing using the M13F and M13R primers. Sequences were visually trimmed to exclude the plasmid backbone and ensure only high-quality sequence data was analyzed. Forward and reverse sequences from individual 470 bp and 1.3 kb insert clones were assembled using the Assemble Sequences tool in CLC Genomics Workbench. Assembled sequences were analyzed using BLASTn to retrieve the most closely related 16S rRNA gene sequences in the NCBI nr database [35].
The Map Reads to References tool in CLC Genomics Workbench was used to map 16S rRNA gene sequences from the Illumina sequencing effort for tanks T30 and T34 (water and sand filter samples) to the Planctomycete-targeted nested PCR sequences. Mapping parameters were set to record sequences with at least 97% sequence identity and a minimum coverage of 50% of the nested PCR sequences. A maximum likelihood phylogenetic tree based on the HKY85 nucleotide substitution model of 16S rRNA gene sequences was generated using PhyML 3.0 software package (http://www.atgc-montpellier.fr/phyml).

Results
As a first step in assessing the presence, diversity, relative abundance, and temporal stability of potential anammox organisms within commercial aquaria, we performed a 16S rRNA gene survey of the microbial communities present in both surface waters and filtration systems for four exhibits within the Tennessee Aquarium. Surface water samples represent the planktonic microbial communities of the aquaria, which are well-mixed and oxygenated. Sand from the filtration systems represent sediment-like environments, within which zones of anoxia are dominant just below the surface. These exhibits, sampled weekly for four consecutive weeks, are characterized by distinct salinities (two marine, two freshwater), temperature profiles (two temperate [~25˚F], two cold [~10˚F and 17˚F]), resident macrofauna, volumes, and turn-over rates (Table 1). Water chemistry measurements were collected in tandem with microbial community sampling to quantitatively compare any fluctuations (S2 Table).

Microbial diversity was stable over time and distinct for each exhibit
A broad-based microbial diversity approach revealed the temporal stability and diversity of microbial communities within these systems. A total of 95 samples were analyzed, yielding 9,912,916 reads passing quality control, resulting in an average of 104,346 reads per sample. A total of 3,845 OTUs (clustered according to >97% identity to cultured organisms) were identified amongst all the samples. Overall, microbial community diversity, as assessed by Shannon-Weiner and Simpson indices, was stable for each habitat over the month-long sampling (Fig  1). With the exception of T34, Shannon-Wiener diversity (H') was significantly higher for communities isolated from sand filters relative to water (Fig 1A). Marine water communities exhibited higher H' values than their freshwater counterparts (Fig 1A). Diversity calculated using the Simpson index (1-λ) also indicated higher diversity in sand filter versus water samples ( Fig 1B). The average Simpson's index of diversity for sand filter samples was 0.95, whereas the average for water samples was 0.77.
Non-metric multidimensional scaling (NMDS) analysis for the one-month sampling period showed that microbial communities clustered according to habitat salinity and substrate-type (Fig 2). Habitats with differing salinity values (marine vs. freshwater) shared �40% similarity (Bray-Curtis metric) in their microbial community compositions. Communities present in the same substrate-type (water vs. sand filter) shared >60% similarity (Bray-Curtis metric). Unique to the marine water samples, temperature was also deterministic of community composition, with communities sharing �69% similarity when grouped by temperature (cold vs. temperate).

Linear discriminant analysis revealed abundant biomarker OTUs
Each habitat had unique dominant microbial OTUs as determined by linear discriminant analysis (LDA) (Fig 3). Sample biomarker OTUs determined by LDA were identified based on 16S rRNA gene abundances that characterized the differences between habitats [36]. For marine samples, the cumulative relative abundance of all biomarkers describing an individual habitat ranged from 13-73% of library reads (Fig 3A and 3B). Relative abundance of freshwater biomarkers ranged from 2-16% of library reads (Fig 3C and 3D). Higher relative abundances of biomarker OTUs were observed in sand filter samples as compared to water samples, regardless of salinity or temperature. In communities derived from the marine temperate water samples (designated T30W), a single OTU (OTU0002), identified as a member of the genus Erythrobacter, dominated the microbial community, comprising 43% of the reads (Fig 3A; S3 Table). Seven biomarker OTUs belonged to the order Planctomycete and were diagnostic for either marine (OTU0008, OTU0025, OTU0029, OTU0185, OTU0295) or freshwater (OTU0010 and OTU0090) sand filter habitats (Fig 4; S3 Table). One sequence type (OTU0008) comprised~10% of the sequence reads from both the cold and temperate marine sand filtration microbial communities (T30S and T34). A second sequence type, OTU0025, comprised~2.5% of these same communities (Fig 3A and 3B; S3 Table).

Metagenomic analyses indicated genetic potential for complete metabolic pathways
Samples of water and sand filtration systems for each of the four habitats, as well as a denitrification system processing water from T30, were used to create metagenomic libraries. From the 9 samples, a total of 30,220,660 reads were obtained. The average contig size for these libraries ranged from~1800-7500 bp (S6 Table), indicating these are diverse microbial communities and a greater depth of sequencing is needed. Despite the relatively low coverage, the metagenome data did indicate the genetic potential for complete pathways of assimilatory sulfate reduction, dissimilatory sulfate reduction and oxidation, and sulfur-oxidation (S2A Fig). The genetic potential for complete pathways of dissimilatory nitrate reduction, assimilatory nitrate reduction, denitrification, and nitrification were also found (S2B Fig). Regarding the anammox pathway, the precursor genes nirK and nirS, which are linked to the conversion of Yet, neither of the essential diagnostic genes for anammox, hydrazine synthase (hzs) or hydrazine dehydrogenase (hdh), were detected [37,38]. Additionally, none of the 16S rRNA genes from reference cultured Planctomycetes isolates mapped to contigs. However, as anammox Planctomycetes that are functionally active often comprise <1% of their community [23,24] and these libraries represented relatively low coverage of these microbial communities, it was reasoned that additional, targeted approaches were needed to more fully address the question of their presence and potential functionality.

Targeted 16S rRNA gene analysis identified multiple Planctomycetes
All samples were further analyzed via Planctomycete-targeted 16S rRNA gene analysis using a nested PCR approach, where the first round of PCR is intended to amplify all Planctomycete sequences and the second round is specific for the anammox-specific sublineage (S1 Fig). However, only sand filter samples from the two marine tanks (T30S and T34S) tanks yielded PCR products. These products were cloned and representatives further analyzed. An initial size survey of 50 clones revealed that <20% of the clones had the anticipated insert size of 500 bp insert; the remaining clones had insert sizes of~1.3 kb, which is equivalent to the product anticipated from the first round of amplification. Ten unique sequences were obtained and the closest matches to sequences in the NCBI database ranged from 88-97% sequence identity (S5 Table). The sequences derived from the longer amplification products had their closest matches to 16S rRNA genes from non-anammox Planctomycetes and analysis of these sequences revealed mismatches with the anammox-specific PCR primer sets employed in the second round of amplification (e.g., S3 Fig). Of the clones representing the size expected for anammox-lineage specific amplification, a single clone (T34-CFU-05) had homology to 16S rRNA genes, with its closest match (88% identity) to a clone denoted as Candidatus 'Brocadia fulgida' (KU217660; S5 Table). Several clones were not 16S rRNA genes, indicative of non-specific amplification. A phylogenetic tree of the clone-derived 16S rRNA sequences supports the homology searches and places T34-CFU-05 near the anammox subgroup within the Planctomycetes (Fig 5). The remaining aquarium-derived sequences from each of the two marine tanks are found on separate branches of the tree.

Discussion
Commercial, large-scale aquaria are engineered to achieve high water clarity, efficient removal of toxic chemical species, and stable environmental conditions to support the resident macrofaunal species. The contributions of resident microorganisms to support a healthy aquatic environment, principally through the transformation of toxic nitrogen species via nitrification and denitrification, has been recognized for over a hundred years [39,40]. However, anammox, a process increasingly employed in wastewater treatment, represents an alternative microbial pathway for processing the most common toxic nitrogen species in aquaria (i.e., ammonium and nitrite), and one that has not been widely considered in this context [41,42]. Anammox Planctomycetes are frequently present in wastewater treatment plants, freshwater lakes, marine suboxic zones, and coastal sediments [20]. The low abundance (typically <1%) of anammox bacteria in natural systems initially prevented identification of these microorganisms. However, the broader understanding of their significant contributions to nitrogen cycling is now well understood and evidence to date suggests group members are present in most, if not all, aquatic systems [20,[22][23][24]. In order to better understand the potential of anammox as a viable process for N control in commercial aquaria, we applied a series of cultureindependent approaches to understand the microbial community structure and potential metabolic functions within Tennessee Aquarium exhibits that differ in chemical and physical properties.
Our broad-based microbial community survey identified Planctomycetes phylum representatives in all aquaria samples. Seven Planctomycete OTUs were identified as diagnostic biomarkers for the different habitats and their relative abundances varied across the microbial communities sampled. These organisms have the greatest representation in both the temperate and cold marine sand filtration systems, where they represent >10% of the microbial communities. Planctomycetes known to carry out anammox appear to be restricted to specific lineages within the phylum and can be best identified by examining variable regions 2 and 3 (V2/V3) of the 16S rRNA gene [43]. Our microbiome analysis utilized the Earth Microbiome primers, which target the V4 region of the 16S rRNA gene. Thus, this approach is not immediately diagnostic of anammox-lineage Planctomycetes. However, it does provide insight into the relative abundance of members of this phylum as well as the full composition and stability of the microbial communities in these systems, all of which are valuable for cross-system comparisons. These types of comparisons will ultimately prove useful in elucidating the role that aquarium microbiomes play in maintaining healthy water chemistry and macrofaunal communities.
In assessing the Tennessee Aquarium metagenomic datasets for genes diagnostic of anammox (i.e., hdh and hzs), neither were identified by sequence homology to functionally confirmed genes. However, two potential precursor proteins (encoded by nirK and nirS) were identified. These proteins encode nitrite reductases that produce ammonia, a key substrate for anammox, but also denitrification [40]. These genes were previously identified in the

PLOS ONE
metagenomes derived from the Georgia Aquarium's Ocean Voyager tank (S2 Fig [17 , 37, 38]). The Georgia Aquarium study assessed the microbial communities residing with sulfur-driven denitrification (SDN) reactors with a focus on anaerobic denitrification performed by sulfuroxidizing bacteria as a method to offload toxic nitrogen species [17]. The SDN systems contained highly diverse microbial communities, and steps of the SDN pathway were predicted to be partitioned amongst community members [17]. All genes involved in nitrogen and sulfur metabolisms identified in Georgia Aquarium Ocean Voyager communities were also found in the Tennessee Aquarium marine metagenomes. Furthermore, the Tennessee Aquarium samples contained additional evidence for nitrogen and sulfur metabolism not evident in the Ocean Voyager metagenomes (S2 Fig; [17]). As was observed in the Ocean Voyager tank [17], the suite of genes required for individual pathways may be distributed across several organisms in the Tennessee Aquarium communities. A more in-depth analysis of the genomic content and phylogenetic context of functional genes is necessary to robustly address this question.
Considering both the absence of functional anammox diagnostic genes and the relatively low sequence coverage of our metagenomes, a 16S rRNA PCR amplification approach specifically targeting anammox Planctomycetes was ultimately employed. This approach provided additional evidence for the presence of Planctomycete phylum representatives within the marine (T30 and T34) sand filter samples (Fig 5). Whether these sequences represent organisms capable of anammox remains unclear. None of the cloned sequences shared >97% identity to any sequences in public databases, and many shared <95% identity (S6 Table). One clone sequence (T34-CFU-5) showed little homology to sequences in the NCBI database, with greatest homology to a candidate anammox lineage member at 88% identity (Fig 5; S4 Table). This low sequence identity is in itself intriguing and further work isolating strains or performing genome assemblies will be needed to characterize this potentially novel bacterium.
The lack of conclusive data supporting the presence of anammox planctomycetes within the Tennessee Aquarium tanks does not necessarily indicate that these bacteria are absent in these environments, nor does it rule out the possibility of bioaugmentation or biostimulation to promote the anammox reaction in these systems. The presence of complete pathways for denitrification and nitrification indicates potential for these habitats to sustain levels of nitrite and ammonium that would enable the anammox reaction [22]. Additionally, wastewater treatment facilities represent systems that contain high concentrations of nitrate, nitrite, and ammonium, yet anammox bacteria are not typically naturally present [44,45]. Instead, they are added to these systems, typically within a bioreactor that provides sufficiently low reduction potential to support the growth of these bacteria [46][47][48][49]. Indeed, Tal et al. (2006), reported successful anammox activity in fixed film biofilters within marine aquaculture systems [49]. These biofilters are structurally and functionally similar to the anammox bioreactors used in wastewater treatment plants [45][46][47][48][49] and represent a path forward for application of analogous systems in commercial aquaria.
Supporting information S1  Fig. Diagram of expected product sizes from nested PCR. The first round of PCR selectively amplifies 16S rRNA genes within the phylum Planctomycete [32,33]. The second round of PCR only amplifies those Planctomycete 16S rRNA genes that belong to the anammox subgroup [34].